
# t3 = in text table 3
# ta8 = online appendix table 8
# fa2 = online appendix fig 2
# ta10 = online appendix table 10 
# ta13 = online appendix table 13 


css2<- read.csv("css05_cv.csv")
css2$part<- css2$cash + css2$vouchers
css2$part[css2$part==2]<-1
css2$part<- css2$part*100
css2$republican<- 0
css2$republican[css2$party=="Democratic"]<- 0
css2$republican[css2$party=="Republican"]<- 1
css2$republican[css2$party=="Non-Partisan"]<- 0
css2$iswhiteperson<- 0
css2$iswhiteperson[css2$race=="European"]<- 1
css2$iswhiteperson[is.na(css2$race)]<- NA
css2$income<- as.numeric(paste0(substr(css2$income,2,nchar(css2$income))))
css2$cash_part<- css2$cash_amt + css2$vouchers_amt


census_blocks2<- css2 %>%
  group_by(bg_ct, elec_year) %>%
  summarize(Mean = mean(part, na.rm=TRUE),
            cash_part = mean(cash_part, na.rm=T))

cb2<- css2 %>% group_by(bg_ct) %>% summarise(pct_white = mean(pct_white, na.rm=T),
                                             iswhiteperson = mean(iswhiteperson, na.rm=T),
                                             income = mean(income, na.rm=T),
                                             med_inc = mean(med_inc,na.rm=T ),
                                             republican = mean(republican, na.rm=T),
                                             sums = mean(sums, na.rm=T))
census_blocks2<- left_join(census_blocks2, cb2)

census_blocks2$pct_white<- census_blocks2$pct_white*100
census_blocks2$iswhiteperson<- census_blocks2$iswhiteperson*100
census_blocks2$republican<- census_blocks2$republican*100
census_blocks2$med_inc<- census_blocks2$med_inc/1000
census_blocks2$income<- census_blocks2$income/1000
census_blocks2$voucher_yr<-0
census_blocks2$voucher_yr[census_blocks2$elec_year>2016]<-1
census_blocks2$q_white<- ceiling(ecdf(census_blocks2$pct_white)(census_blocks2$pct_white)*4)
census_blocks2$rep_q<- ceiling(ecdf(census_blocks2$republican)(census_blocks2$republican)*4)
census_blocks2$sums_q<- ceiling(ecdf(census_blocks2$sums)(census_blocks2$sums)*4)


census_blocks3<- css2 %>%
  group_by(bg_ct19, elec_year) %>%
  summarize(Mean = mean(part, na.rm=TRUE),
            cash_part = mean(cash_part, na.rm=T))

cb3<- css2 %>% group_by(bg_ct19) %>% summarise(
  income = mean(income, na.rm=T),
  med_inc = mean(med_inc,na.rm=T ))
census_blocks3<- left_join(census_blocks3, cb3)
census_blocks3$med_inc<- census_blocks3$med_inc/1000
census_blocks3$income<- census_blocks3$income/1000
census_blocks3$voucher_yr<-0
census_blocks3$voucher_yr[census_blocks3$elec_year>2016]<-1
census_blocks3$inc_q<- ceiling(ecdf(census_blocks3$med_inc)(census_blocks3$med_inc)*4)

se_white<-summary(lm_robust(Mean~voucher_yr*pct_white, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= census_blocks2)) #t3

se_whiteq<-summary(lm_robust(Mean~voucher_yr*as.factor(q_white), 
                             fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                             se_type="stata", clusters=bg_ct, data= census_blocks2))# ta10
 
#pct republican
se_rep<-summary(lm_robust(Mean~voucher_yr*republican, 
                          fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                          se_type="stata", clusters=bg_ct, data= census_blocks2)) #t3
se_repq<-summary(lm_robust(Mean~voucher_yr*as.factor(rep_q), 
                           fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                           se_type="stata", clusters=bg_ct, data= census_blocks2))# ta10

#l2 pct white
se_iwp<-summary(lm_robust(Mean~voucher_yr*iswhiteperson, 
                          fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                          se_type="stata", clusters=bg_ct, data= census_blocks2)) #ta8

lm1<- lm(pct_white~iswhiteperson, data=census_blocks2)
ggplotRegression <- function (fit) {
  require(ggplot2)
  ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
    geom_point() +
    stat_smooth(method = "lm", col = "red") +
    labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
                       " Slope =",signif(fit$coef[[2]], 5),
                       " P =",signif(summary(fit)$coef[2,4], 5)))+
    xlab("L2 pct white")+ylab("Census block pct white")
}
white_vs_l2<-ggplotRegression(lm1) # fa2


#income


se_blockinc<-summary(lm_robust(Mean~voucher_yr*med_inc, 
                               fixed_effects=~as.factor(elec_year)+as.factor(bg_ct19), 
                               se_type="stata", clusters=bg_ct19, data= census_blocks3)) #t3
se_blockincq<-summary(lm_robust(Mean~voucher_yr*as.factor(inc_q), 
                                fixed_effects=~as.factor(elec_year)+as.factor(bg_ct19), 
                                se_type="stata", clusters=bg_ct19, data= census_blocks3))# ta10

#income with l2
se_l2inc<-summary(lm_robust(Mean~voucher_yr*income, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct19), 
                            se_type="stata", clusters=bg_ct19, data= census_blocks3)) #ta8

lm1<- lm(med_inc~income, data=census_blocks3)
ggplotRegression <- function (fit) {
  
  require(ggplot2)
  
  ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
    geom_point() +
    stat_smooth(method = "lm", col = "red") +
    labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
                       " Slope =",signif(fit$coef[[2]], 5),
                       " P =",signif(summary(fit)$coef[2,4], 5)))+
    xlab("L2 Income")+ylab("Census block median income")
}
inc_vs_l2<-ggplotRegression(lm1) # fa2


#mean voting record by tract

se_bvote<-summary(lm_robust(Mean~voucher_yr*sums, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= census_blocks2)) #t3

se_bvoteq<-summary(lm_robust(Mean~voucher_yr*as.factor(sums_q), 
                             fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                             se_type="stata", clusters=bg_ct, data= census_blocks2))# ta10

#with cash now 
se_white_c<-summary(lm_robust(cash_part~voucher_yr*pct_white, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= census_blocks2)) #ta13
se_bvote_c<-summary(lm_robust(cash_part~voucher_yr*sums, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= census_blocks2)) #ta13
se_blockinc_c<-summary(lm_robust(cash_part~voucher_yr*med_inc, 
                               fixed_effects=~as.factor(elec_year)+as.factor(bg_ct19), 
                               se_type="stata", clusters=bg_ct19, data= census_blocks3)) #ta13
se_rep_c<-summary(lm_robust(cash_part~voucher_yr*republican, 
                          fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                          se_type="stata", clusters=bg_ct, data= census_blocks2)) #ta13



#placebo test (table a27)
cb2<- subset(census_blocks2, census_blocks2$elec_year<2016)
cb2$voucher_yr[cb2$elec_year>2012]<-1
cb3<- subset(census_blocks3, census_blocks3$elec_year<2016)
cb3$voucher_yr[cb3$elec_year>2012]<-1

se_white_pl<-summary(lm_robust(Mean~voucher_yr*pct_white, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= cb2)) 

se_rep<-summary(lm_robust(Mean~voucher_yr*republican, 
                          fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                          se_type="stata", clusters=bg_ct, data= cb2)) 
se_blockinc<-summary(lm_robust(Mean~voucher_yr*med_inc, 
                               fixed_effects=~as.factor(elec_year)+as.factor(bg_ct19), 
                               se_type="stata", clusters=bg_ct19, data= cb3)) 
se_bvote<-summary(lm_robust(Mean~voucher_yr*sums, 
                            fixed_effects=~as.factor(elec_year)+as.factor(bg_ct), 
                            se_type="stata", clusters=bg_ct, data= cb2)) 
